home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / dgbmv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  11.2 KB  |  361 lines
  1. *
  2. ************************************************************************
  3. *
  4. *     File of the DOUBLE PRECISION  Level-2 BLAS.
  5. *     ===========================================
  6. *
  7. *     SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
  8. *    $                   BETA, Y, INCY )
  9. *
  10. *     SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
  11. *    $                   BETA, Y, INCY )
  12. *
  13. *     SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
  14. *    $                   BETA, Y, INCY )
  15. *
  16. *     SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
  17. *    $                   BETA, Y, INCY )
  18. *
  19. *     SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
  20. *
  21. *     SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  22. *
  23. *     SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
  24. *
  25. *     SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
  26. *
  27. *     SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  28. *
  29. *     SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
  30. *
  31. *     SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
  32. *
  33. *     SUBROUTINE DGER  ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  34. *
  35. *     SUBROUTINE DSYR  ( UPLO, N, ALPHA, X, INCX, A, LDA )
  36. *
  37. *     SUBROUTINE DSPR  ( UPLO, N, ALPHA, X, INCX, AP )
  38. *
  39. *     SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  40. *
  41. *     SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
  42. *
  43. *     See:
  44. *
  45. *        Dongarra J. J., Du Croz J. J., Hammarling S.  and Hanson R. J..
  46. *        An  extended  set of Fortran  Basic Linear Algebra Subprograms.
  47. *
  48. *        Technical  Memoranda  Nos. 41 (revision 3) and 81,  Mathematics
  49. *        and  Computer Science  Division,  Argonne  National Laboratory,
  50. *        9700 South Cass Avenue, Argonne, Illinois 60439, US.
  51. *
  52. *        Or
  53. *
  54. *        NAG  Technical Reports TR3/87 and TR4/87,  Numerical Algorithms
  55. *        Group  Ltd.,  NAG  Central  Office,  256  Banbury  Road, Oxford
  56. *        OX2 7DE, UK,  and  Numerical Algorithms Group Inc.,  1101  31st
  57. *        Street,  Suite 100,  Downers Grove,  Illinois 60515-1263,  USA.
  58. *
  59. ************************************************************************
  60. *
  61.       SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
  62.      $                   BETA, Y, INCY )
  63. *     .. Scalar Arguments ..
  64.       DOUBLE PRECISION   ALPHA, BETA
  65.       INTEGER            INCX, INCY, KL, KU, LDA, M, N
  66.       CHARACTER*1        TRANS
  67. *     .. Array Arguments ..
  68.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  69. *     ..
  70. *
  71. *  Purpose
  72. *  =======
  73. *
  74. *  DGBMV  performs one of the matrix-vector operations
  75. *
  76. *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
  77. *
  78. *  where alpha and beta are scalars, x and y are vectors and A is an
  79. *  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
  80. *
  81. *  Parameters
  82. *  ==========
  83. *
  84. *  TRANS  - CHARACTER*1.
  85. *           On entry, TRANS specifies the operation to be performed as
  86. *           follows:
  87. *
  88. *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  89. *
  90. *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  91. *
  92. *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
  93. *
  94. *           Unchanged on exit.
  95. *
  96. *  M      - INTEGER.
  97. *           On entry, M specifies the number of rows of the matrix A.
  98. *           M must be at least zero.
  99. *           Unchanged on exit.
  100. *
  101. *  N      - INTEGER.
  102. *           On entry, N specifies the number of columns of the matrix A.
  103. *           N must be at least zero.
  104. *           Unchanged on exit.
  105. *
  106. *  KL     - INTEGER.
  107. *           On entry, KL specifies the number of sub-diagonals of the
  108. *           matrix A. KL must satisfy  0 .le. KL.
  109. *           Unchanged on exit.
  110. *
  111. *  KU     - INTEGER.
  112. *           On entry, KU specifies the number of super-diagonals of the
  113. *           matrix A. KU must satisfy  0 .le. KU.
  114. *           Unchanged on exit.
  115. *
  116. *  ALPHA  - DOUBLE PRECISION.
  117. *           On entry, ALPHA specifies the scalar alpha.
  118. *           Unchanged on exit.
  119. *
  120. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  121. *           Before entry, the leading ( kl + ku + 1 ) by n part of the
  122. *           array A must contain the matrix of coefficients, supplied
  123. *           column by column, with the leading diagonal of the matrix in
  124. *           row ( ku + 1 ) of the array, the first super-diagonal
  125. *           starting at position 2 in row ku, the first sub-diagonal
  126. *           starting at position 1 in row ( ku + 2 ), and so on.
  127. *           Elements in the array A that do not correspond to elements
  128. *           in the band matrix (such as the top left ku by ku triangle)
  129. *           are not referenced.
  130. *           The following program segment will transfer a band matrix
  131. *           from conventional full matrix storage to band storage:
  132. *
  133. *                 DO 20, J = 1, N
  134. *                    K = KU + 1 - J
  135. *                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
  136. *                       A( K + I, J ) = matrix( I, J )
  137. *              10    CONTINUE
  138. *              20 CONTINUE
  139. *
  140. *           Unchanged on exit.
  141. *
  142. *  LDA    - INTEGER.
  143. *           On entry, LDA specifies the first dimension of A as declared
  144. *           in the calling (sub) program. LDA must be at least
  145. *           ( kl + ku + 1 ).
  146. *           Unchanged on exit.
  147. *
  148. *  X      - DOUBLE PRECISION array of DIMENSION at least
  149. *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  150. *           and at least
  151. *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  152. *           Before entry, the incremented array X must contain the
  153. *           vector x.
  154. *           Unchanged on exit.
  155. *
  156. *  INCX   - INTEGER.
  157. *           On entry, INCX specifies the increment for the elements of
  158. *           X. INCX must not be zero.
  159. *           Unchanged on exit.
  160. *
  161. *  BETA   - DOUBLE PRECISION.
  162. *           On entry, BETA specifies the scalar beta. When BETA is
  163. *           supplied as zero then Y need not be set on input.
  164. *           Unchanged on exit.
  165. *
  166. *  Y      - DOUBLE PRECISION array of DIMENSION at least
  167. *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  168. *           and at least
  169. *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  170. *           Before entry, the incremented array Y must contain the
  171. *           vector y. On exit, Y is overwritten by the updated vector y.
  172. *
  173. *  INCY   - INTEGER.
  174. *           On entry, INCY specifies the increment for the elements of
  175. *           Y. INCY must not be zero.
  176. *           Unchanged on exit.
  177. *
  178. *
  179. *  Level 2 Blas routine.
  180. *
  181. *  -- Written on 22-October-1986.
  182. *     Jack Dongarra, Argonne National Lab.
  183. *     Jeremy Du Croz, Nag Central Office.
  184. *     Sven Hammarling, Nag Central Office.
  185. *     Richard Hanson, Sandia National Labs.
  186. *
  187. *     .. Parameters ..
  188.       DOUBLE PRECISION   ONE         , ZERO
  189.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  190. *     .. Local Scalars ..
  191.       DOUBLE PRECISION   TEMP
  192.       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY,
  193.      $                   LENX, LENY
  194. *     .. External Functions ..
  195.       LOGICAL            LSAME
  196.       EXTERNAL           LSAME
  197. *     .. External Subroutines ..
  198.       EXTERNAL           XERBLA
  199. *     .. Intrinsic Functions ..
  200.       INTRINSIC          MAX, MIN
  201. *     ..
  202. *     .. Executable Statements ..
  203. *
  204. *     Test the input parameters.
  205. *
  206.       INFO = 0
  207.       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
  208.      $         .NOT.LSAME( TRANS, 'T' ).AND.
  209.      $         .NOT.LSAME( TRANS, 'C' )      )THEN
  210.          INFO = 1
  211.       ELSE IF( M.LT.0 )THEN
  212.          INFO = 2
  213.       ELSE IF( N.LT.0 )THEN
  214.          INFO = 3
  215.       ELSE IF( KL.LT.0 )THEN
  216.          INFO = 4
  217.       ELSE IF( KU.LT.0 )THEN
  218.          INFO = 5
  219.       ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN
  220.          INFO = 8
  221.       ELSE IF( INCX.EQ.0 )THEN
  222.          INFO = 10
  223.       ELSE IF( INCY.EQ.0 )THEN
  224.          INFO = 13
  225.       END IF
  226.       IF( INFO.NE.0 )THEN
  227.          CALL XERBLA( 'DGBMV ', INFO )
  228.          RETURN
  229.       END IF
  230. *
  231. *     Quick return if possible.
  232. *
  233.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  234.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  235.      $   RETURN
  236. *
  237. *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  238. *     up the start points in  X  and  Y.
  239. *
  240.       IF( LSAME( TRANS, 'N' ) )THEN
  241.          LENX = N
  242.          LENY = M
  243.       ELSE
  244.          LENX = M
  245.          LENY = N
  246.       END IF
  247.       IF( INCX.GT.0 )THEN
  248.          KX = 1
  249.       ELSE
  250.          KX = 1 - ( LENX - 1 )*INCX
  251.       END IF
  252.       IF( INCY.GT.0 )THEN
  253.          KY = 1
  254.       ELSE
  255.          KY = 1 - ( LENY - 1 )*INCY
  256.       END IF
  257. *
  258. *     Start the operations. In this version the elements of A are
  259. *     accessed sequentially with one pass through the band part of A.
  260. *
  261. *     First form  y := beta*y.
  262. *
  263.       IF( BETA.NE.ONE )THEN
  264.          IF( INCY.EQ.1 )THEN
  265.             IF( BETA.EQ.ZERO )THEN
  266.                DO 10, I = 1, LENY
  267.                   Y( I ) = ZERO
  268.    10          CONTINUE
  269.             ELSE
  270.                DO 20, I = 1, LENY
  271.                   Y( I ) = BETA*Y( I )
  272.    20          CONTINUE
  273.             END IF
  274.          ELSE
  275.             IY = KY
  276.             IF( BETA.EQ.ZERO )THEN
  277.                DO 30, I = 1, LENY
  278.                   Y( IY ) = ZERO
  279.                   IY      = IY   + INCY
  280.    30          CONTINUE
  281.             ELSE
  282.                DO 40, I = 1, LENY
  283.                   Y( IY ) = BETA*Y( IY )
  284.                   IY      = IY           + INCY
  285.    40          CONTINUE
  286.             END IF
  287.          END IF
  288.       END IF
  289.       IF( ALPHA.EQ.ZERO )
  290.      $   RETURN
  291.       KUP1 = KU + 1
  292.       IF( LSAME( TRANS, 'N' ) )THEN
  293. *
  294. *        Form  y := alpha*A*x + y.
  295. *
  296.          JX = KX
  297.          IF( INCY.EQ.1 )THEN
  298.             DO 60, J = 1, N
  299.                IF( X( JX ).NE.ZERO )THEN
  300.                   TEMP = ALPHA*X( JX )
  301.                   K    = KUP1 - J
  302.                   DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL )
  303.                      Y( I ) = Y( I ) + TEMP*A( K + I, J )
  304.    50             CONTINUE
  305.                END IF
  306.                JX = JX + INCX
  307.    60       CONTINUE
  308.          ELSE
  309.             DO 80, J = 1, N
  310.                IF( X( JX ).NE.ZERO )THEN
  311.                   TEMP = ALPHA*X( JX )
  312.                   IY   = KY
  313.                   K    = KUP1 - J
  314.                   DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL )
  315.                      Y( IY ) = Y( IY ) + TEMP*A( K + I, J )
  316.                      IY      = IY      + INCY
  317.    70             CONTINUE
  318.                END IF
  319.                JX = JX + INCX
  320.                IF( J.GT.KU )
  321.      $            KY = KY + INCY
  322.    80       CONTINUE
  323.          END IF
  324.       ELSE
  325. *
  326. *        Form  y := alpha*A'*x + y.
  327. *
  328.          JY = KY
  329.          IF( INCX.EQ.1 )THEN
  330.             DO 100, J = 1, N
  331.                TEMP = ZERO
  332.                K    = KUP1 - J
  333.                DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL )
  334.                   TEMP = TEMP + A( K + I, J )*X( I )
  335.    90          CONTINUE
  336.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  337.                JY      = JY      + INCY
  338.   100       CONTINUE
  339.          ELSE
  340.             DO 120, J = 1, N
  341.                TEMP = ZERO
  342.                IX   = KX
  343.                K    = KUP1 - J
  344.                DO 110, I = MAX( 1, J - KU ), MIN( M, J + KL )
  345.                   TEMP = TEMP + A( K + I, J )*X( IX )
  346.                   IX   = IX   + INCX
  347.   110          CONTINUE
  348.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  349.                JY      = JY      + INCY
  350.                IF( J.GT.KU )
  351.      $            KX = KX + INCX
  352.   120       CONTINUE
  353.          END IF
  354.       END IF
  355. *
  356.       RETURN
  357. *
  358. *     End of DGBMV .
  359. *
  360.       END
  361.